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Strong evidence for the Helling-Downs correlations have been reported by sev- 
eral pulsar timing array collaborations in middle 2023. In this work, we study the 
state-of-the-art graviton mass bounds by analyzing the observational data of overlap 
ruduction functions from NANOGrav 15-year data release and CPTA first data re- 
lease. The data analysis places upper limits on the graviton mass at 95% confidence 
level, namely, mg < 0.43 x 107*3eV for NANOGrav and mg < 0.57 x 1077%eV for 


CPTA. In addition, we discuss implications of these results for scenarios of ultralight 


tensor dark matter. 


I. INTRODUCTION 


Based on Einstein’s theory of general relativity, the Hellings-Downs (HD) correlation 
curves have been proposed to characterize a stochastic gravitational-wave background 
(SGWB) of pulsar timing array (PTA) band [1]. Recently, strong evidence for a stochastic 
signal that is spatially correlated among multiple pulsars were reported by the North Amer- 
ican Nanohertz Observatory for Gravitational Waves (NANOGrav) [2] and Chinese PTA 
(CPTA) [3] Collaborations, respectively. In particular, the HD correlations were claimed 


to deserve statistical significance of ~ 30 — 40 by NANOGrav and of 4.60 by CPTA. The 
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European PTA (EPTA) [4] and Parkes PTA (PPTA) [5] Collaborations also claimed that 
their datasets are compatible with the HD correlations. 

However, the overlap reduction functions (ORFs) for a theory of massive gravity could 
be different from the HD correlations. They depend on the graviton speed and then on the 
graviton mass, but recovers the HD correlations in the massless limit [6, 7]. The massive 
gravity was first proposed by Fierz and Pauli in 1939 [8]. It has been extensively studied in 
the subsequent more than eighty years and the graviton mass has been bounded by a quantity 
of laboratory and astronomy observations (e.g., see reviews in Ref. [9] and references therein). 

During the era of gravitational-wave astronomy, the speed of gravitational waves v, was 
shown by the Advanced LIGO, Virgo and KAGRA (LVK) Collaborations to almost coincide 
with the light speed with a precision |v, — 1| < 107'°, indicating an upper limit on the 
graviton mass, i.e., my < 107?” eV [10]. Recently, similar upper bounds were also claimed 
by other research groups |11, 12], who focused on analysis of the past NANOGrav 12.5-year 
dataset [13]. 

In this work, taking into account ORFs for massive gravitons in Ref. [7], we explore 
the state-of-the-art upper bounds on the graviton mass by analyzing the 2023 NANOGrav 
and CPTA datasets |2, 3] that supported strong evidence for the HD correlations but did 
not veto other alternative correlations. The remaining context of the paper is as follows. In 
Section I, we show ORFs for massive gravitons, with the HD correlations being the massless 
limit. In Section II, we show a likelihood method for data analysis as well as the resulting 


bounds on the graviton mass. In Section IV, we show the concluding remarks. 


I. THEORY 


In this work, we can consider only the helicity-two modes for simplicity [14]. Conven- 
tionally, we may study the massive gravity by analyzing its Lagrangian. However, when we 
study the propagation only, an equivalent but simple approach seems to concern the disper- 
sion relation, i.e., W? = k? + m2, with the graviton speed being defined as vy = dw/dk [10]. 
Here, m, is the graviton mass, w and k stand for the angular frequency and wavenumber, 
respectively. For subluminal gravitons, we obtain a deviation of the graviton speed from the 


light speed, i.e., 


w= 1- (™) (1) 


W 


For a typical frequency band, the lower limit on the graviton speed can be recast into the 
upper bound on the graviton mass. Throughout this work, the light speed is defined as 
unity. 

Upon a SGWB, the time of arrivals of radio pulses from two pulsars would be spatially 


correlated, with the angular correlation defined as 


Yab(Vg) = 5 t LoP Cab) ; (2) 
€ 


where the subscript ap stands for the cross correlation of two pulsars a and b with the angular 


separation Cab, and the angular power spectrum Cy is defined as 
Ce = T"? Je(vg, f Da) JE (vg f Do) , (3) 


where f = w/(27) is the frequency of gravitational waves, De denotes a distance to the c-th 


pulsar. To simplify the above definition, we have introduced a new function of the form 


o | (EQ)! [7s dx om, Jela) 
Keay) = Vmi ay fe (4) 


where je(x) denotes the spherical Bessel function with -th multipole. Following Refs. [7, 15], 


we recast the angular correlation 7,,(v,) into ORFs P',,(v,) by normalizing the former such 
that Paa(v, = 1) = 0.5, where the subscript aa stands for the auto correlation of a-th pulsar. 
Note that the HD correlation curves are recovered by the above correlations in the massless 
limit. 

In Fig. | and Fig. 2, we show the difference of ORFs between the massive gravitons and the 
massless gravitons. For comparison, we also reproduce the observed ORFs by NANOGrav 
[2] (see Fig. 1) and CPTA [3] (see Fig. 2). Roughly, it seems that the observed ORFs fit 
the massive gravitons better than the massless gravitons. The data analysis in the following 


section would be proved to support such suspicions. 


II. DATA ANALYSIS AND RESULTS 


For the NANOGrav 15-year data release, there are 67 pulsars with a timing baseline longer 
than 3 years monitored, implying 2,211 distinct pairs in total. Each pair has a deterministic 
angular separation. Based on those datasets, the authors of Ref. [2] constructed a minimally 


modeled Bayesian reconstruction of the inter-pulsar correlation pattern by employing spline 
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FIG. 1. Comparison between the ORFs of massive gravitons and those of massless gravitons. The 


orange shaded region stands for the cosmic variance. The NANOGrav 15-year data points are 


shown as shaded violins for comparison. 
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FIG. 2. Comparison between the ORFs of massive gravitons and those of massless gravitons. The 


orange shaded region stands for th cosmic variance. The CPTA first-release data points are shown 


for comparison. 


interpolation over seven spline-knot positions, i.e., Fig. 1(d) of the paper. For the CPTA first 


data release [3], there are 57 millisecond pulsars monitored. The 4.60 statistical significance 


of the HD correlation between those pulsars has been found around 14 nHz, see Fig. 4 of 


the data release paper. 


Analyzing the above NANOGrav and CPTA datasets of the spatial correlations, we infer 
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FIG. 3. Posteriors of the graviton speed vg inferred from the NANOGrav 15-year data release and 


the CPTA first data release. The vertical dashed lines stand for the 95% confident lower limits. 


the graviton speed, or equivalently, the graviton mass, by using the log-likelihood as follows 


—2In L(v,[D) = > (220 T) i (5) 


Cab Cab 


where TP, stands for the observed ORFs for the binned angular separation Cab, 7, stands 
for the corresponding lo uncertainty, D denotes the NANOGrav or CPTA datasets, and the 
summation runs over all the binned angular separations. Note that the above uncertainties 
also contain the cosmic variance [7, 15, 16]. 

For our data analysis, the parameter to be inferred is the graviton speed vg. We let the 
corresponding priors to be uniform, i.e., vg € [10~?, 1], and conduct a Markov-Chain Monte- 
Carlo sampling by using the public code cobaya [17]. We also use the public PTAfast package 
[15] to compute ORFs. The resulting posteriors of v, will be recast into the posteriors of 
mg, following the relation in Eq. (1). Therefore, we will then obtain the upper limit on mg 
at 95% confidence level. 

Our results are shown as follows. The posteriors of vg are depicted in Fig. 3. We label the 
95% confident lower limits on vg with vertical dashed lines. They are shown to be vg Z 0.13 
for NANOGrav and vy Z 0.14 for CPTA. Compared with NANOGrayv, CPTA favors a 


~~ 


relatively smaller speed of gravitons. We further display the posteriors of m, in Fig. 4. The 
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FIG. 4. Posteriors of the graviton mass inferred from the NANOGrav 15-year data release and 
the CPTA first data release. The shaded regions stand for the allowed parameter space at 95% 


confidence level. 


shaded regions stand for the allowed parameter space at 95% confidence level. Therefore, we 
find upper limits on m, at 95% confidence level, namely, mg < 4.3x 107?3 eV for NANOGrav 
and m, < 5.7 x 107” eV for CPTA. We find that the two bounds are consistent with each 


other. They stand for the state-of-the-art upper bounds on the graviton mass. 


IV. CONCLUSION 


In this work, we investigated the graviton mass bounds through analysis of the NANOGrav 
15-year dataset and the CPTA first data release. Analyzing the data points of ORFs ob- 
served by the two observatories, we inferred the allowed parameter interval for the graviton 
speed, particularly, the posteriors. Recasting the posteriors of graviton speed into the pos- 
teriors of graviton mass, we obtained the state-of-the-art upper limits on the graviton mass, 
i.e., Mg S few x 10-78 eV at 95% confidence level. 

There may be implications of our results for the ultralight tensor dark matter [18, 19], 


which may account for the mystery of dark matter in the universe. In particular, the 


ultralight dark matter in the mass range Mulidm ~ 10~22eV was used for conquering several 
shortcomings of the traditional cold dark matter [20, 21]. The ultralight tensor dark matter 
behaves like the massive gravitons, indicating possible imprints of it on PTAs. Since our 
results were obtained through analysis of the recent PTA datasets, we find that the preferred 
mass range of ultralight tensor dark matter is compatible with the upper bounds on graviton 
mass of this current work. A similar work resulting from analysis of the NANOGrav 12.5- 


year dataset can be found in Ref. [22]. 
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